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ABSTRACT 

Reliable measurements of the solar magnetic field are still restricted to the photo- 
sphere, and our present knowledge of the three-dimensional coronal magnetic field is 
largely based on extrapolation from photospheric magnetogram using physical mod- 
els, e.g., the nonlinear force-free field (NLFFF) model as usually adopted. Most of 
the currently available NLFFF codes have been developed with computational vol- 
ume like Cartesian box or spherical wedge while a global full-sphere extrapolation is 
still under developing. A high-performance global extrapolation code is in particu- 
lar urgently needed considering that Solar Dynamics Observatory (SDO) can provide 
full-disk magnetogram with resolution up to 4096 x 4096. In this work, we present a 
new parallelized code for global NLFFF extrapolation with the photosphere magne- 
togram as input. The method is based on magnetohydrodynamics relaxation approach, 
the CESE-MHD numerical scheme and a Yin- Yang spherical grid that is used to over- 
come the polar problems of the standard spherical grid. The code is validated by 
two full-sphere force-free solutions from Low & Lou's semi-analytic force-free field 
model. The code shows high accuracy and fast convergence, and can be ready for 
future practical application if combined with an adaptive mesh refinement technique. 

Subject headings: Magnetic fields; Magnetohydrodynamics (MHD); Methods: nu- 
merical; Sun: corona 



1. Introduction 

Magnetic field holds a central position within solar research such as sunspots and coronal 
loops, prominences, solar flares, solar wind and coronal mass ejections. However a routinely di- 
rect measurement of solar magnetic field that we can rely on is restricted to the solar surface, i.e., 
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the photosphere, in spite of the works that have been done to measure the coronal fields using the 
radio and infrared wave bands ( Gary & Hurford|1994| Lin et al.|2004 ). This is extremely unfortu- 
nate since the magnetic field plays a comparatively minor role in the photosphere but completely 
dominates proceedings in the corona ( |Solanki et al. 2006| ). Up to present, our knowledge of the 
three-dimensional (3D) coronal magnetic field is largely based on extrapolations from photo spheric 
magnetograms using physical models. For the low corona, one model assuming free of Lorentz 
force (J x B = 0, where J is the current and B is the magnetic field) is justified by a rather small 
plasma (3 (the ratio of gas pressure to magnetic pressure) and a quasi-static state. The force-free 
assumption involves a intrinsically nonlinear equations (V x B) x B = that is rather difficult to be 
solved if based on boundary information alone, and various computing codes have been proposed 
to solve this equation numerically for nonlinear force-free field (NLFFF) extrapolations (e.g., see 
review papers by |Amari et al] ( |1997| ), |McClymont et al.| ( |1997| ), |Schrijver et al.| ( |2006[ ), |Metcalf 
etaL] ( 12555) , |Wiegelmann| ( [2008] ) and |DeRosa etaT|P^ ). 



Most of the currently available NLFFF codes are developed in Cartesian coordinates. Thus 
the extrapolations are limited to relatively local and small areas, e.g., a single active region (AR) 
without any relationship with other ARs. However, the ARs usually cannot be isolated since they 
interact with the neighboring ARs or overlying large-scale fields. Observations of moving plasma 
connecting several separated ARs by SOHO Extreme Ultraviolet Imaging Telescope (EIT) reveal 
the connections between ARs (e.g., [Wang et al. ( 2001| )). Also the activities in the chromosphere 
and corona often spread over several ARs, such as filament bursts recorded in Ha images and 
coronal mass ejections (CMEs) observed by SOHO Large Angle and Spectrometric Coronograph 
(LASCO) coronagraphs. Even for a single AR, it is pointed out that the fields of view in Cartesian 



box are often too small to properly characterize the entire relevant current system (DeRosa et al. 



2009). To study the connectivity between multi-ARs and extrapolate in a larger field of view, it 
is necessary to take into account the curvature of the Sun's surface by extrapolation in spherical 
geometry partly or even entirely, i.e., including the global corona ( Wiegelma nn|2007" Tadesse et al 



2011 2012[ ). Moreover, a global NLFFF extrapolation can avoid any lateral artificial boundaries 
which cause issues in Cartesian codes. Global non-potential extrapolation are urgently needed 
considering that high-resolution, full-disk vector magnetograms will be soon available from the 
Solar Dynamics Observatory (SDO). Another motivation comes from the developing of global 
MHD models for the solar corona and solar wind. Up to present the global MHD models are based 
on only the line-of-sight (LoS) magnetogram, using the global potential extrapolation to initialize 



the computation (e.g., Feng et al. (2010)). These models will be challenged by the full-sphere 
vector magnetogram, which need a global non-potential extrapolation. 

In the past few years, several global NLFFF extrapolation methods have been developed but 
they are still in its infancy and many issues need to be resolved. For example: |He & Wang| (2006) 
validate the boundary-integral-equation method (Yan & Li 2006| ) for extrapolation above a full 
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sphere using simple models of |Low & Lou| ( |1990j ) while application to more complex extrap- 
olations needs further development; Wiegelmann (2007) and Tadesse et al.| ( [2009] ) extend their 
optimization code to spherical coordinates including for both partial and full sphere, but the con- 



vergence speed is proved to be rather slow if polar regions are included in the computation ( Wiegel- 
mann|[2007 ); A flux rope insertion method based on magnetofriction has also been developed by 



van B allego oijen| ([2004) for constructing NLFFFs in spherical coordinates (e.g., see its applica- 
tions by |Bobra et al.| ( [2008] ); |Su et al.| ( |2009a|b[ ); |Savcheva & van Ballegoorjen] ( [2009] )), but the 
same problem as in Wie gelmann] ( |2007| ) maybe encounted if the code is extended to containing the 
whole sphere; very recently [Contopoulos et al. (2011 ) present a new force-free electrodynamics 
method for global coronal field extrapolation, however the solution is not unique since it is only 
prescribed by the radial magnetogram. 

Among the existing problems, choosing a suitable grid system is in particular critical for the 
implementation of any global models (not only the global force-free extrapolation) with the lower 
boundary as the full solar surface. Naturally one can use the spherical grid, i.e., with grid lines 
defined by coordinates (r, 9,(f)). However, the simplicity of a standard- spherical coordinate grid is 
destroyed by the problem of grid convergence and grid singularity at both poles (Usma nov]|1996 



Kageya ma & Sato||2004[ [Feng et aT7||2010[ ). These problems severely restrict the converge speed 



of the full-sphere optimization code as reported by Wiegelmann (2007). Although the singular- 
ity problem can be partially resolved by excluding a small high-latitude cone, e.g., restricting the 
latitude within 10° < 9 < 170°, it inescapably disconnects the field lines crossing over the polar 
region. To avoid these problems, Contopoulos et al. ( 20 1 1 [ ) simply use the Cartesian grid to imple- 
ment their method by taking a cubic Cartesian box to contain the whole region with the solar sphere 
cut out. But the Cartesian grid can not characterize precisely the Sun's sphere surface, since the 
solar sphere nowhere coincides with any grid points. Thus the corresponding boundary conditions 
is hard to prescribe. On the other hand, the unstructured grid has been used frequently in global 



MHD models ( |Tanaka|[1994| |Feng et aL]|2007| |Nakamizo et al.|[2009| ) and can possibly be intro- 
duced into NLFFF extrapolation, but it costs heavily in mesh generation and management and is 
not suitable for solvers based on numerical difference (but may be suitable for finite element solver, 



e.g., the FEMQ code developed by Amari et al. (2006)). Moreover on the unstructured grid, it is 
difficult to implement the technique of parallelized-adaptive mesh refinement (AMR), which is an 
attractive tool for resolving the contradiction between the computational demand for extrapolating 
high-resolution/large field-of-view magnetograms and the computational resource limitations. A 
promising solution to the above problems is use of an overlapping spherical grid, e.g., see several 
types of overlapping spherical grid proposed by |Usmanov| ( |1996| ), |Kageyama & Sato| ( |2004[ ),[Hen- 
shaw & Schwendeman| ( |2008] ) and |Feng et al.| ( |2010] ). In principle one can use a set of low-latitude 
partial-sphere grids to cover the full sphere with some patches overlapped. Certainly in such over- 
lapping grid system, there is no pole problems and meanwhile the grid management is easy. If the 
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component grids are carefully chosen, the overlapping patches can be minimized and only add a 
very small numerical overhead on the computation for data communications between the compo- 
nent grids. Among the overlapping grids for composing the full sphere, a so-called Yin- Yang grid 
(Kageya ma & Sato|2004| ) is a most elegant configuration with only two identical components and 



the overlapping region less than 7% of the full sphere (see Figure[T|). 



Our previous work <\ Jiang et aT|2011| | Jiang & Feng||2012j ) has been devoted to a new imple 



mentation of MHD relaxation approach for NLFFF extrapolation based on the CESE numerical 
scheme (space-time conservation-element and solution-element scheme). We have introduced a 
new set of magneto-frictional-like equations with the AMR and a multigrid-like strategy for accel- 
erating the computation and improving its convergence. The good performance and high accuracy 
of the code has been demonstrated by detailed comparisons with previous work by |Schrijver et al. 



( |2006[ ) and |Metcalf et al.| ((2008 ) based on several NLFFF benchmark tests, although it remains 



to be seen how well this method will work with real solar data. The success of the CESE-MHD- 
NLFFF code encourages us to extend it to spherical geometry and ultimately realize a fast and 
accurate way for global non-potential extrapolation of the high-resolution full-disk magnetograms 
from SDO. In this paper, we take the first step by developing a new code for global NLFFF extrap- 
olation using the CESE-MHD-NLFFF method on the Yin- Yang grid. This code is assumed to be 
applied to force-free vector magnetograms, i.e., with the uncertainties and inconsistencies removed 



by some kind of preprocessing approach, e.g., that proposed by Wiegelmann et al. (2006). This 
is important considering that NLFFF codes generally failed to extrapolate a satisfactory force-free 
field when applied to non-force-free magnetogram (Met calf et al.||2008| ), which is usually the case 
of observed data. 

The remainder of the paper is organized as follows. In Section |2]we give the model equations 
and the numerical method including a curvilinear-version CESE-MHD scheme and its implemen- 
tation in Yin- Yang grid. In Section [3] we set up two semi- analytic test cases of full-sphere NLFFF 
solution proposed by Low & Lou| ( |1990 ). The extrapolation results and qualitative and quantitative 



comparisons are presented in Section |4} Finally, we draw conclusions and give some outlooks for 
future work in Section |5] 



2. The Method 

2.1. Model Equations 

The basic idea of using the MHD relaxation approach to solve the force-free field is to use 
some kind of fictitious dissipation to drive the MHD system to an equilibrium in which all the 
forces can be neglected if compared to the Lorentz force and the boundary vector map is satisfied. 
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In our previous work for NLFFF extrapolation (Jiang & Feng 2012), a magnetic splitting form of 



magneto-frictional model equations is introduced as 

^ = (VxB 1 )xB-z/pv, 

at 

<9B 

— - = Vx(vxB)+V(uV'B 1 )-vV-B 1) 
at 

art 

-^ = 0,VxBo = 0,V-Bo = 0, 

at 

p=|B| 2 +p , 6 = 60+6!. (1) 

In the above equation system: 6 is a potential field matching the normal component of the mag- 
netogram, 6! is the deviation between the potential field 6 and the force-free field 6 to be solved; 
v is the frictional coefficient and p is a numerical diffusive speed of the magnetic monopole; po is 
a necessary small value (e.g., p = 0.01) to deal with very weak field associated with the magnetic 
null. The value for parameters v and p are respectively given by v = At/ Ax 2 and p = 0.4 Ax 2 / At, 
according to the time step At and local grid size Ax. Merits of using the above equations include: 

• Retaining explicitly the time-dependent form of the momentum equatiorQmake the magneto- 
frictional equation system able to be handled by modern CFD or MHD solver designed for 
the standard partial-differential-equation system like 

dU | d¥(V) | <9G(U) | <9H(U) =s 

dt dx dy dz 

• Numerically, accuracy can be gained for solving only the deviation field 61 by dividing the 
total magnetic field 6 into two parts (6 = 6 + 6i) (Tanaka [1994 ). Also such a splitting has a 



physical meaning ( Priest & Forbes|2002[ ): potential component 6 arises from photospheric 



or sub-photospheric currents and can be regarded as invariant during a flare, whereas the non- 
potential component B\ arises from large-scale coronal currents (above the photosphere) and 
is the source of the flare energy. 

Any numerical magnetic monopoles V ■ Bi can be rapidly convected away with the plasma 
by term -vV • 61 and effectively diffused out by term V(pV ■ 61). 

Setting a pseudo-plasma density p oc |6| 2 can equalize the Alfven speed of the whole domain 
and thus accelerate the relaxation in the weak field regions. 



This is unlike the standard magneto-frictional method in which the momentum equation is simplified as vy = J x B. 
Although the standard magneto-frictional method is simpler since only the induction equation ^ = V x (v x B) = 
V x ( 1^5 x B) is needed to be solved, this equation cannot be written in the form of Equation (IZb. 
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2.2. Numerical Implementation 



To solve the model equation ([T]) in spherical geometry, we employ a curvilinear version of the 
CESE-MHD solver proposed by Jiang et al. ( 2010[ ). In this method, the governing equations written 
as Equation ([2]) are transformed from the physical space (x, y, z) to a reference space (£, 77, Q whose 
mapping is explicitly known x = x(£, 77, Q;j =;y(£, rj, Q;z = z(£, rj, Q. The transformed equations are 



where 



dU OF dG <9H 

dt <9£ drj dC, 



U = 7U, 

F=/(Fe,+Ge y +He z ), 

G = 7(Fr/ ;c + Gr /v + Hr/ z ), 
H = 7(FC, + GC>. + HQ, 
S = 7S; 



(3) 



(4) 



and 7 is determinant of Jacobian matrix J for the mapping, i.e., 



d(x,y,z) 



y$ y v yc 

Z£ Z v Z{ 



(5) 



Based on this transformation, the basic idea is to map the spherical geometry of physical space 
to a simple rectangular grid of the reference space, in which the we can use the Cartesian CESE- 
MHD method to solve the transformed equations with very simple rectangular-uniform mesh. For 



detailed descriptions of the CESE-MHD method and its curvilinear version please refer to (Jiang 
et al.|2010|[Feng et al.|2012l ). 



To overcome the grid- singularity problem at both poles of the standard spherical coordinates, 
we use the Yin- Yang grid. As a type of overlapping grid, the Yin- Yang grid is synthesized by 
two identical component grids in a complemental way to cover an entire spherical surface with 
partial overlap on their boundaries (see Figure [T]). Each component grid is a low latitude part of 
the latitude-longitude grid without the pole. Therefore the grid spacing on the sphere surface is 
quasi-uniform and the metric tensors (i.e., the matrix elements in Equation ([5])) are simple and 
analytically known (Kageya ma & Sato||2004| ). In Figure [T] one component grid, say 'Yin' grid, is 
defined in the spherical coordinates by 



\6-ir/2\ <7r/4+5;|0-7r| <3n/4+5 



(6) 



where 5 = 1.5A0 is a small buffer to minimize the required overlap. The other component grid, 
'Yang' grid, is defined by the same rule of Equation ([6]) but in another coordinate system that is 
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rotated from the Yin's, and the relation between Yin coordinates and Yang coordinates is denoted 
in Cartesian coordinates of each own by (x e ,y e ,z e ) = (-x n ,z n ,y„), where (x n ,y n ,z n ) is Yin's Cartesian 
coordinates and (x ei y e ,z e ) is Yang's. 

We then map the Yin and Yang component grids to rectangular grids by defining two mapping 
equations 

{x = e^ sin#cos0 
y = e^s\vi9sm(j) (7) 
z = e^ cos6> 

and 

{x = -e^s'm9cos(p 
y = e^cos9 , (8) 

z = sin^sin^ 

where (£, 9, <p) is the coordinates of the reference space with rectangular-uniform mesh used ( A£ = 
A9 = A0). In this definition, we have 

Ar = e 5+Af - <K = e 5 (e A? rA£ = rA# (9) 

which means that the cells are close to regular cubes in physical space, especially at low latitudes. 

The grid extent in £ is £ e [0,ln 10], i.e., the outer boundary is set at r = lORs (solar radius). 
The initial condition is specified by simply setting B\ = and v = 0. The constant part B is 
obtained by a fast potential field solver which is developed by a combination of the spectral and 
the finite-difference methods (Jiang & Feng 2012| ). The lower boundary condition is given by 



the vector magnetogram while the outer boundary is fixed with zero values of B t and v. In the 
following test cases, we focus the field extrapolation in r e [l,2]i? s and the objective of setting the 
outer boundary far beyond the extrapolation volume is to minimize the boundary effect. 

On the boundaries where grids overlap, solution values on one component grid are deter- 
mined by interpolation from the other. We use explicit interpolation for simplicity and effi- 
ciency in parallel computation, and the grid buffer 5 is suitably chosen for enough overlap area 
to perform such interpolation (see Figure [T]). In the reference space, standard tensor-product 
Lagrange interpolation ( jlsaacson & Keller 1966j ) is used. For instance (see Figure [2] for de- 



tails), the interpolation of values / at the point M(£ m ,Vm,(m) in the reference space is computed 
by KM) = £Lo£?=oE?=o/(U>fcyf (OP? (v)Pk (0. where P?(x) is the Lagrange interpolating 
polynomial P?(x) = Ylk=o k=/j 73^ w i m x being £, 77 or (. Note that the interpolation accuracy is of 
three order which is higher than the CESE solver by one order. Thus the discretization accuracy in 
the overlapping region will not be reduced by the interpolation. Finally, to realize the paralleliza- 
tion on this bi-component grids, each component grid is divided into small blocks, e.g., consisting 
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of 8 x 8 x 8 cells with guard-cells (one layer of ghost cells for convenience of communication 
between blocks), which are distributed evenly among the processors. Message-passing-interface 
(MPI) library is employed for data communications between the processors. The interpolation of 
the overlapping boundaries is dealt with in a similar way as for the intra-grid guard-cell filling and 
both operations are arranged to be done simultaneously. The load balancing is also considered 
carefully among all the processors to further improve the parallel scaling. 



3. Test Case 



The NLFFF model derived by |Low & Lou| ([T990) has served as standard benchmark for many 
extrapolation codes ( |Wheatland et al.|2000[ |Amari et al.l|2006t |Schrijver et~aLl|2006t |Valori et al.' 
2007 1 He & Wang|2008 Jiang et al.|2011| ). The fields of this model are basically axially symmetric 
and can be represented by a second-order ordinary differential equation of P(/S) derived in spherical 
coordinates 



d 2 P 



l + n 



(\-H 2 ) — + n(n+\)P + a 2 P l+1/n = 

dfi 1 n 



(10) 



where n and a are constants and /i = cos#. With boundary conditions of P = at fx = — 1, 1, the 



solution P of Equation ( 10) is uniquely determined by two eigenvalues, n and its number of nodes 
m ( Low & Lou|1990 ; Amari et al.|2006 ). The magnetic fields are then given by 



B r = 



1 dA 
r 2 sin# 09 



1 dA 
rsinO dr 



Brk — 



1 



rsin^ 



Q 



(ll) 



where A = P(fi)/r n and Q = aA l+l l n . The fields are axisymmetric in spherical coordinates (i.e., 
invariant in <p direction) with a point source at the origin. To avoid such obvious symmetry in 
the M1-3D extrapolation test, we locate the point source with 03R S offset to the center of the 
computational volume and deviate the axis of symmetry with the z-axis by $ = 7r/10 (the length 
unit of the above equations is the solar radius R s ). We present two test cases with eigenvalues 
of n = l,m = 1 (hereafter referred to as CASE LL1) and n = 3,m = 1 (CASE LL2) respectively. 
Both cases are performed on the same resolution of 90 x 180 grids in 9-(j) plane and r e [l,2]R s . 
The synoptic maps of the field and the force-free parameter a at the bottom of the solutions are 
shown in Figures [3] and [4] and the 3D field lines are shown in panel (a) of Figures [5] and |6j It 
should be remarked that these magnetograms do not represent the real magnetic distributions of the 
photosphere but only be used for the purpose of testing our code. As can be seen, the a distribution 
of CASE LL2 is more inhomogeneous than that of CASE LL1, which means that CASE LL2 is 



more nonlinear. We note that CASE LL1 is very similar to test cases used by Wiegelmann (2007 1 



and Tad esse et al.| p009) while CASE LL2 is more difficult than tests in their works. 



Before inputting the vector maps in the NLFFF code, we made some consistency checks for 
the maps. If a vector map is used for a force-free extrapolation, some necessary conditions have to 
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be fulfilled ( |Aly||1989 Sakurai 1989, Tadesse et al. 2009). For clarification we repeat here these 
conditions from Tadesse (201 1 1 where a detailed derivation of the condition formula is given. First 
of all the net magnetic flux must be in balance, i.e., 



B r ds = 



(12) 



where S represents the whole sphere. Secondly the total force on the boundary has to vanish, which 
can be expressed in spherical coordinates as 



Ty = 



1 



00 +B 2 b -B 2 ) sin 9 cos - B r Bg cos 9 cos + B r B^ sin < 



1 



(Bg + BI-B 2 ) sin 9 sin - B r B e cos 9 sin - B,B,p cos < 



■F 3 = 



s L 



1 



{B 2 e +B 2 ) -B 2 r )cos9 + B r B e sin9 



ds = 0; 
ds = 0; 
ds = 0. 



(13) 



Thirdly the total torque on the boundary vanishes, i.e., 



7i = / BjiB^ cos 9 cos <p+B e sin 0) ds = 0; 
is 

?2 = B r (B<p cos 9 sin <f)-B e cos <p) ds = 0; 
Js 



T 3 = jB r B,sm9ds = Q. 



(14) 



To quantify the quality of the synthetic full-disk magnetograms with respect to the above criteria, 
we compute three parameters, i.e., the flux balance parameter 



Eflu 



the force balance parameter 
and the torque balance parameter 



S S \Br\ ds' 
[J 7 !] + 1^21 + |^3 1 

E B 

\Ti\ + \T 2 \ + 1751 



orque 



(15) 



(16) 



(17) 



where E B = J^B^+Bl+B^) ds. For the above cases, the three parameters are (-7.1 17 x 10 4 , 1 . 181 x 
10" 4 , 6.313 x 10" 6 ) and (-1.028 x 10" 3 , 5.458 x 10" 4 , 7.379 x 10" 5 ) respectively, which shows that 
these maps are ideally consistent with the force-free model. 
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4. Results 

In this section, we present the results of extrapolation and compare them with their original 
solutions qualitatively and quantitatively. As usual, the quantitative comparison is performed by 
computing a suite of metrics (also referred to as figures of merit), which are listed as follows: 

• the vector correlation C vec 

Cve^^Brb^^lB,-! 2 ^^-! 2 ), (18) 



• the Cauchy-Schwarz inequality C C s 



• the normalized and mean vector error E n , E m 



F = — V" l B '~ b 'l . J7 1 _ 1 _ F 

^ m ~M^ IB-I ' m_ 

• the magnetic energy ratio e 



(21) 



V lb-1 2 

e= W <22) 

where B ; and b, denote the Low & Lou solution and the extrapolated field, respectively, i denotes 
the indices of the grid points and M is the total number of grid points involved. It is also important 
to measure the ratio of the total energy to the potential energy 

<23) 

to study the free energy budget for realistic coronal field. 

We also calculate another four metrics to measure the force-freeness and divergence-freeness 
of the results. They are the current- weighted sine metric CWsin 

CWsin^S^;^^^, (24 ) 
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the divergence metric (\fi\) 



and the £ V . B and £ V xB 



j_ v (V-B),. 
KUlW Af^6|B,|/Ax' 



(25) 



_ _ 1 v-. iBiCV-B),-] 



Af-VlV(|B| 2 /2),r 



Ui X B; 



M^|V(|B|V2),. 



(26) 



All the above metrics have been described in detail in our previous work Jiang & Feng (2012) and 
thus will not be repeated here. 



4.1. Qualitative Comparison 

In Figures[5]and|6} we present a side-by-side comparisons of the extrapolation results with the 
Low & Lou models and the potential fields by plotting the 3D field configurations. For each cases 
the field lines are traced from the same set of footpoints on the photosphere. A good agreement 
between our extrapolation results and the original Low & Lou solutions can be seen from the 
highly similarity of most of the field lines. The basic difference between the extrapolated fields 
and the potential fields is the shearing, which is reconstructed by the bottom-boundary-driving 
process exerted on the initial un-sheared potential fields. By placing the outer boundary far away 
enough, we can make most of the field lines move freely in the volume, which thus is helpful for 
the relaxation of the field lines. Figure [7] compares the field values of the r = 1.57? 5 surface by 
plotting the contours of the reference solutions (solid lines) and the extrapolation (dashed lines) on 
the same figure. As can be seen, contours lines of the fields from the reference solution and the 
extrapolation almost overlapped with each other. 



4.2. Quantitative Comparison 

Quantitative metrics shown in Table [T] and [2] demonstrate good performance of the code. In 
these tables, we present results of the full sphere with r E [l,2]R s and more lower region r E 
[1, 1.5]i? s . For both cases results of the vector correlation C vec and Cqs are extremely close to the 
reference values, showing a perfect matching of the vector direction. Results of vector error E' n and 
E' m also score close to 1 (even the error of most sensitive metric E' m is smaller than 10%), while 
the potential solutions have results of only ~ 0.5. This is very encouraging, since in previously 
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reported tests of Cartesian or spherical NLFFF extrapolation code with only photospheric boundary 
provided, e.g., done by |Schrijver et al.| ( |2006[ ); [Valori et al.] ( |2007[ ); [Wiegel mann| ( |2007[ ); |Tadesse 



et al. (2009), results with E' m > 0.9 are rarely achieved. These two metrics show that the original 



solutions are reconstructed with very high accuracy. Finally, the energy content of the non-potential 
fields, a critical parameter from the extrapolation used to calculate the energy budget in solar 
eruptions, is also well reproduced (with errors under several percents). By comparison of the 
metrics, we find that accuracy of the lower region is even higher than the full region, which means 
the strong fields are extrapolated better than the upper-region weak fields. In real solar fields, only 
the lower part of the corona is close to force-free while the upper corona is no more force-free 
because of the expansion of hot plasma ( |Gary|[200T ). Thus to extrapolate the lower-region fields 
more close to the original force-free solution than the upper-region fields is consistent with the real 
solar conditions. 

Table [3] gives the metrics measuring force-freeness and divergence-freeness of the fields, 
which are rather small and close to the level of discretization error. Unlike the first two metrics 
(CWsin and (|//|)) which mainly characterize the geometric properties of the field, metrics £VxB 
and £ V b are introduced to measure the physical effect of the residual divergence and Lorentz 



forces on the system in the actual numerical computation (Jiang & Feng 2012). This is impor- 
tant when checking the NLFFF solution if it is used to initiate any full-MHD simulations. Of the 
present extrapolation results, the residual forces are less than one percent of the magnetic-pressure 
force. 



4.3. Convergence Study 



We finally give a study of convergence of the extrapolation. In Figure [8] and [9] we show how 
the system relaxes and reaches its final force-free equilibrium by plotting the temporal evolution 
of several parameters, including the residual of field Bi 



res"(Bi) = 



\ 



=x,y,z 



£«) 2 



(27) 



(where n denotes the iteration steps^J, the maximum and average velocity, and the nine metrics 
described above. The system converged very fast from a initial residual of > 10" 1 to value ~ 10" 5 
with time of IOOta (about 6000 ~ 8000 iteration steps, see panel (a) of the figures). The evolution 



2 in the present experiments, we record the parameters not every time step but every eight steps for saving the total 
computing time, thus n-8 means the residual is for the full eight steps. 
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of the plasma velocity indicates that initially (1) the system is driven away from the starting v = 
state and then (2) by the relaxation process a static equilibrium is reached as expected with a rather 
small residual velocity which is only on the order of the numerical error 0(Ax 2 ) of the CESE 
solver. All the metrics plotted in the figures converged after A0r A (less than 2000 iteration steps), 
when the residual is on the order of 10" 4 , and the convergence speed of CASE LL2 is even faster 
than CASE LL1. Note that the metrics (\fj\) and £y. B , like the plasma velocity, first climbs to a 
relatively high level (see panel (d)) and then drops to the level of discretization errors. In principle 
the divergence-free constraint of B should be fulfilled throughout the evolution, at least close to 
level of discretization error. However, an ideally dissipationless induction equation Equation ([T]) 
with divergence-free constraint can preserve the magnetic connectivity, which makes the topology 



of the magnetic field unchangeable (Wiegelmann 2008) unless a finite resistivity is included to 



allow the reconnection and changing of the magnetic topology ( Roumeliotis|| 1 996 ) . In the present 
implementation in which no resistivity is included in the induction equation, an allowing of high 
values in V ■ B in the initial evolution process (indicated by the climb of metric (\fi\)) may provide 
some freedom for changes in the magnetic topology (also note that a numerical diffusion can help 
topology adjustment). 

Besides the extrapolation accuracy, the computing time also matters for a practical use of the 
NLFFF code. In this present tests, the size of the Yin- Yang grid is equivalent to 80 x 90 x 180 in 
ordinary spherical grid. The computation is completed by less than 2 hours using 32 processors 
of Intel Xeon CPU E5450 (3.00GHz). In practical applications, the computing time can be further 
reduced considering that it is not necessary to evolve the system to 100r A . 



5. Conclusions 

In this work we present a new code for NLFFF extrapolation of the global corona. The method 
is implemented by installing the previous code CESE-MHD-NLFFF in the Cartesian geometry 
onto a Yin- Yang spherical grid. By this grid system, we can incorporate intrinsically the full- 
sphere computation and avoid totally the problems involved with the spherical poles. The boundary 
conditions are only specified on the bottom sphere and free of any lateral-boundary information. 
We have examined the performance of this newly developed code using two test cases of the 
classic semi-analytic force-free fields by Low & Lou ( 1990). We show that the code runs fast and 



achieves a good accuracy with the extrapolation solution very close to the reference field and the 
force-freeness and divergence-freeness constraints well fulfilled. 

We note that the success of extrapolating the model solutions (i.e., the ideal force-free test 
cases) does not necessarily indicate successful applications to the real solar data, which contains 
various inconsistencies and uncertainties. A practicable solution to this issue may be attributed 
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to some preprocessing methods as developed by Wiegelmann et al. (2006) and Fuhrmann et al 



( |2007 ) or more realistic MHD model focused on the photosphere-chromosphere interface with at- 



mosphere stratification. For another important issue, reconstruction that need to be performed with 
much larger field of view (to more entirely characterize the currents between ARs and high over 
the AR) and with higher spatial resolution (to capture the fine critical structures such magnetic 
null point) is severely limited by the computational capability. Application with grid of about 500 3 
pixels is almost the upper limit for computational capability of most recently developed or updated 
codes. This is rather unsatisfied if considering extrapolation with the 4096 x 4096 SDO/HMI mag- 
netograms. This issue can be resolved promisingly via combining the global extrapolation code 
with the adaptive mesh refinement, according to the intrinsic characteristics of the solar magnetic 
field in which the active regions represents only a small fraction of the whole surface. By the 
AMR technique, one can take focus on local corona, e.g., some active regions in the context of 
global extrapolation with the corresponding high-resolution vector magnetograms embedded in a 
low-resolution global vector or LoS map, as exemplified by Figure 10 In this figure f\ mesh for 



the active regions are refined with three more grid levels than the background full-sphere grid (a 
block- AMR algorithm (Powell et al. 1999[ ) is used and the ratio of resolutions between the grid 



levels is two). The grid structure can be dynamically adjusted during the MHD-relaxation process 
to capture the strong currents and important magnetic structures (e.g., flux ropes) by designing 
carefully the AMR refinement criteria. Our future work will include performing more stringent 
testing of the code and installation of the code on AMR grid for practical application to SDO/HMI 
data. 
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Model 


C V ec 


Ccs 


F' 


F' 


e 


E/E po t 


For r e [1,2]7? 5 
Low 

Extrapolation 
Potential 


1 

0.9995 
0.8595 


1 

0.9974 
0.8204 


1 

0.9609 
0.5261 


1 

0.9269 
0.4641 


1 

0.9783 
0.8517 


1.1741 
1.1486 
1 


For re [l,1.5]fl s 
Low 

Extrapolation 
Potential 


1 

0.9998 
0.8620 


1 

0.9995 
0.8236 


1 

0.9772 
0.5441 


1 

0.9668 
0.5013 


1 

0.9851 
0.8780 


1.1390 
1.1220 
1 


Table 1: CASE LL1: Results of the metrics. 


Model 


C V ec 


Ccs 


F' 


F' 


e 


E/Ep t 


Forr e [l,2]fl s 
Low 

Extrapolation 
Potential 


1 

0.9999 
0.9049 


1 

0.9965 
0.8013 


1 

0.9807 
0.5733 


1 

0.9456 
0.4515 


1 

1.0061 
0.9056 


1.1042 
1.1110 
1 


For re [l,1.5]fl s 
Low 

Extrapolation 
Potential 


1 

0.9999 
0.9055 


1 

0.9997 
0.8529 


1 

0.9864 
0.5861 


1 

0.9810 
0.5174 


1 

1.0063 
0.9092 


1.0999 
1.1068 
1 



Table 2: CASE LL2: Results of the metrics. 



Case 


CWsin 


<l//|> 


£VxB 




LL1 


2.11 x 10" 2 


3.73 x 10" 4 


9.04 x 10" 3 


1.83 x 10" 2 


LL2 


2.68 x 10" 2 


4.09 x 10" 4 


5.94 x 10" 3 


7.64 x 10" 3 



Table 3: Results for the metrics measuring the force-freeness and divergence-freeness. 
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Fig. 2. — Interpolation of the point M in the reference space: M represents a mesh point on the 
overlapping boundary, for example, of Yin grid, then the twenty-seven (3 3 ) nodes denote the Yang's 
inner mesh points that are closest to M. 
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Fig. 3. — Vector map (upper) and a distribution (bottom) of CASE LL1. In the vector map, 
contours represent B r and the tangential field (B^,Bq) is shown by the vectors with blue color 
positive B r region and red in negative B r region. 
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Fig. 4.— Same as Figure|3]but for CASE LL2. 
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Fig. 5. — CASE LL1: 3D magnetic field lines with contour of B r on the photosphere surface, (a) 
Low & Lou's solution; (b) the extrapolation result; (c) the initial potential field. 




Fig. 6.— Same as Figure|5]but for CASE LL2. 
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Fig. 7. — Contour map comparison of the Low & Lou's solution (solid lines) and extrapolation 
solution (dashed lines) at r = l.5Rs for B r (top), Bg (middle) and (bottom). Left column is for 
CASE LL1 and right for CASE LL2. 
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Fig. 8. — CASE LL1: The history of the relaxation to force-free equilibrium, (a) Evolution of 
residual resCBO with time (and the iteration steps); (b) evolution of the maximum and average 
velocity; and (c), (d) evolution of the metrics. 



-26- 




Fig. 10. — A example of global extrapolation with AMR grid. The sphere represents the pho- 
tosphere with grid structure outlined by the black lines (note that here each mesh cell represent 
one grid block). The color lines represent the magnetic field lines and the grey image represents 
the radial field B r . A slice of z = -0.3 is plotted to show the radial mesh structure. Note that the 
high-resolution blocks are clustered around the active regions. 



